# -*- coding: utf-8 -*-
"""
Calculation of ballistics of titan lander (with heat-protective coating burning and
without it)
"""

import os
import shutil

import matplotlib.pyplot as plt
import numpy as np
from docx.shared import Cm
from docxtpl import DocxTemplate, InlineImage
from scipy.integrate import ode
from scipy.optimize import fsolve
from titanatmosphere import h_atm, interpolate, rho_atm

G = 6.6743e-20  # Гравитационная постоянная, км^3/(кг*с^2)

r_titan = 2576  # Средний радиус Титана, км
a_titan = 1221870  # Большая полуось орбиты Титана, км
e_titan = 0.0288  # Эксцентриситет орбиты Титана
M_titan = 1.3452 * 10**23  # Масса Титана, кг
f_titan = G * M_titan  # Гравитационный параметр Титана, км^3/c^2
g_titan = 0.001352  # Ускорение свободного падения Титана, км/c^2

spline = interpolate.make_interp_spline(h_atm, rho_atm)
f = f_titan * 10**9
R = r_titan * 10**3

# Параметры опорной орбиты
a0 = R + 2000 * 10**3  # Большая полуось опорной орбиты, м
e0 = 0.2  # Эксцентриситет опорной орбиты
w0 = 0.05  # Аргумент перигея опорной орбиты
p0 = a0 * (1 - e0**2)  # Перицентр опорной орбиты, м

# Параметры точки входа в атмосферу; скорость Титана в перигее - 5.576 км/с
theta_vh = -20 * np.pi / 180  # Угол входа
psi_vh = 0 * np.pi / 180  # Широта точки входа
h_vh = 900 * 10**3  # Высота точки входа, м
vI = np.sqrt(f / (R + h_vh))  # I космическая скорость, м/с
v_vh = vI + 0.0  # Скорость входа, м/с
r_vh = R + h_vh  # Радиус точки входа, м

# Параметры спускаемого аппарата
r_sa = 2.7 / 2  # Радиус СА, м
Sm = np.pi * r_sa**2  # Площадь миделя, м^2
Cx = 1  # Коэффициент А/Д сопротивления
K = 0.3  # Коэффициент А/Д качества
Cy = K * Cx  # Коэффициент подъемной силы
M = 320  # Масса СА, кг

Td_T = 427 + 273
rho_T = 2210
eta_T = (8337.2 + 7117.5) * 10**3 / 2  # все для тефлона

# Характеристики ТЗП
Td_tz = Td_T
rho_zt = rho_T
eta_tz = eta_T


Ck = 2.24 * 10**8
Cr = 0.435
Rn = 1.5 * np.sqrt(Sm / np.pi)
eps = 0.95
const_sb = 5.67e-8

# Диапазон рассматриваемых значений
YSm = [0.50, 0.75, 1, 1.25, 1.5]  # Множитель миделя
Ypsi_vh = [0, 120, 240]  # Перечень широт точки входа, град
Yv_vh = [-300, -150, 0, 150, 300]  # Слагаемое к I космической скорости, м/с
Ytheta_vh = [-75, -60, -45, -30, -15]  # Перечень углов входа, град
YK = [0.10, 0.25, 0.40, 0.55, 0.70, 0.85, 1.00, 1.15]  # Перечень K
YCx = [0.15, 0.2, 0.4, 0.6, 0.8, 1.00, 1.2, 1.4]  # Перечень Cx


def ballistic(v, theta, psi):
    global p1, a1, e1, w1, psi_sh, dV, gamma

    # параметры орбиты
    C1 = r_vh * v * np.cos(theta)
    p1 = C1**2 / f
    H1 = v**2 - 2 * f / r_vh
    a1 = -f / H1
    e1 = np.sqrt(1 - p1 / a1)
    u_vh = np.arccos((p1 - r_vh) / (e1 * r_vh))
    w1 = psi - u_vh

    # Аргументы точки схода
    def func(psi_sh):
        return p0 / (1 + e0 * np.cos(psi_sh - w0)) - p1 / (1 + e1 * np.cos(psi_sh - w1))

    psi_sh = fsolve(func, psi)
    u0_sh = psi_sh - w0
    u1_sh = psi_sh - w1

    # Компоненты скоростей
    Vn0 = np.sqrt(f / a0) * (1 + e0 * np.cos(u0_sh)) / np.sqrt(1 - e0**2)
    Vr0 = np.sqrt(f / a0) * (e0 * np.sin(u0_sh)) / np.sqrt(1 - e0**2)
    Vn1 = np.sqrt(f / a1) * (1 + e1 * np.cos(u1_sh)) / np.sqrt(1 - e1**2)
    Vr1 = np.sqrt(f / a1) * (e1 * np.sin(u1_sh)) / np.sqrt(1 - e1**2)
    print(Vn0, Vr0, Vn1, Vr1)

    # Тормозной импульс и его компоненты, угол приложения
    dVn = abs(Vn1 - Vn0)
    dVr = abs(Vr1 - Vr0)
    dV = np.sqrt(dVn**2 + dVr**2)
    gamma = np.arcsin(dVr / dV)
    print(dVn, dVr, dV, np.degrees(gamma))
    return p1, a1, e1, w1, u_vh, psi_sh, u0_sh, u1_sh, dV, gamma


# Обработчик шага для случая без учета уноса ТЗП
Y1, T1, t_list1, y_list1 = [], [], [], []


def step1(t, y):
    global Y1, T1
    V, theta, L, h = y
    t_list1.append(t)
    y_list1.append(list(y.copy()))
    if h <= 100:
        Y1 = np.array(y_list1)
        T1 = np.array(t_list1)
        return -1


# Обработчик шага для случая с учетом уноса ТЗП
Y2, T2, t_list2, y_list2 = [], [], [], []


def step2(t, y):
    global Y2, T2
    V, theta, L, h, m = y
    t_list2.append(t)
    y_list2.append(list(y.copy()))
    if h <= 10:
        Y2 = np.array(y_list2)
        T2 = np.array(t_list2)
        return -1


# Функция правых частей СДУ для случая без учета уноса ТЗП
def equations1(t, y):
    V, theta, L, h = y
    return [
        (-Cx * Sm * spline(h) * V**2 / 2 - M * f * np.sin(theta) / (h + R) ** 2) / M,
        (
            Cy * Sm * spline(h) * V**2 / 2
            - M * f * np.cos(theta) / (h + R) ** 2
            + M * V**2 * np.cos(theta) / (h + R)
        )
        / (M * V),
        V * np.cos(theta) * R / (h + R),
        V * np.sin(theta),
    ]


# Функция правых частей СДУ для случая с учетом уноса ТЗП
def equations2(t, y):
    V, theta, L, h, m = y
    T_op = Tp(h, V) >= Td_tz
    return [
        -Cx * Sm / (M - T_op * Sm * m) * spline(h) * V**2 / 2
        - f * np.sin(theta) / (h + R) ** 2,
        Cy * Sm / (M - T_op * Sm * m) * spline(h) * V / 2
        - f * np.cos(theta) / (V * (h + R) ** 2)
        + V * np.cos(theta) / (h + R),
        V * np.cos(theta) * R / (h + R),  # В лекциях R/(h+R) нет, но должно быть
        V * np.sin(theta),
        T_op * (qk(h, V) + qr(h, V) - eps * const_sb * Td_tz**4) / eta_tz,
    ]


# Уравнения тепловых потоков и температуры поверхности ТЗП
def qk(h, V):
    return Ck * (spline(h) / spline(0)) ** 0.5 * (V / vI) ** 3.2 / np.sqrt(Rn)


def qr(h, V):
    return Cr * Rn * (spline(h) / spline(0)) ** 1.168 * (V / vI) ** 7.4


def Tp(h, V):
    return ((qk(h, V) + Sm * qr(h, V)) / (eps * const_sb)) ** 0.25


# Определение полярного угла
def Yfi_find(T, fi, Yv, Yh):
    k = 0
    fi_list = [fi]
    for k in range(len(T) - 1):
        tau1, rv1, Vel = T[k], R + Yh[k], Yv[k]
        k = k + 1
        tau2, rv2 = T[k], R + Yh[k]
        tau = tau2 - tau1
        rv = (rv2 + rv1) / 2
        lan = tau * Vel
        fi = fi - lan / rv
        fi_list.append(fi)
    return fi_list


# Определение максимума коэффициента перегрузки
imax = 0


def Yn_nmax(Cx, Sm, Yh, Yv, K):
    global imax
    Yn = (
        Cx
        * Sm
        * spline(Yh)
        * Yv**2
        / (2 * M * f / (Yh + R) ** 2)
        * np.sqrt(1 + K**2)
    )
    i = 0
    imax = 0
    while i <= len(Yh) - 2:
        if Yn[i - 1] < Yn[i] > Yn[i + 1] and Yn[i] > Yn[imax]:
            imax = i
        i = i + 1
    return Yn


# Нахождение максимума температуры и начала/завершения уноса ТЗП
k1, k2, kmax = 0, 0, 0


def temp_points(T, Yh, Yv):
    global k1, k2, kmax
    k = 0
    k_check = 0
    while k <= len(T) - 2:
        if (
            Td_tz > Tp(Yh[k], Yv[k]) < Tp(Yh[k + 1], Yv[k + 1])
            and Tp(Yh[k + 1], Yv[k + 1]) > Td_tz
            and k_check == 0
        ):
            k1 = k + 1
            k_check += 1
        if Td_tz < Tp(Yh[k], Yv[k]) and Tp(Yh[k - 1], Yv[k - 1]) > Td_tz:
            k2 = k
        if Tp(Yh[k - 1], Yv[k - 1]) < Tp(Yh[k], Yv[k]) > Tp(
            Yh[k + 1], Yv[k + 1]
        ) and Tp(Yh[k], Yv[k]) > Tp(Yh[kmax], Yv[kmax]):
            kmax = k
        k = k + 1


# BLOK_VISUALISACII


# Shema_POLETA
def flight_scheme(title, Yh, Yfi):
    """Function risovaniya trajectorii poleta

    Parameters
    ----------
    title : string
        nazvanie grafika
    Yh : list of float
        visota poleta (atmospherniy uchastok)
    Yfi : list of float
        ugol poleta (atmospherniy uchastok)
    """
    fig1, ax1 = plt.subplots()
    fig1.set_facecolor("white")
    plt.axis("equal")
    angle = np.linspace(0, 2 * np.pi, 1000)
    plt.plot(
        p0 / (1 + e0 * np.cos(angle - w0)) * np.cos(angle) * 10 ** (-3),
        p0 / (1 + e0 * np.cos(angle - w0)) * np.sin(angle) * 10 ** (-3),
        color="black",
        linewidth=1.5,
        linestyle="dashed",
        label="Опорная орбита",
    )
    plt.plot(
        (R + h_vh) * np.cos(angle) * 10 ** (-3),
        (R + h_vh) * np.sin(angle) * 10 ** (-3),
        linestyle=":",
        linewidth=0.75,
        color="black",
        label="Граница атмосферы",
    )
    plt.plot(
        p1 / (1 + e1 * np.cos(angle - w1)) * np.cos(angle) * 10 ** (-3),
        p1 / (1 + e1 * np.cos(angle - w1)) * np.sin(angle) * 10 ** (-3),
        color="black",
        linewidth=0.75,
        linestyle="dashed",
        label="Орбита спуска",
    )
    plt.plot(
        (Yh + R) * np.cos(Yfi) * 10 ** (-3),
        (Yh + R) * np.sin(Yfi) * 10 ** (-3),
        linestyle="solid",
        linewidth=1,
        color="black",
        label="Траектория спуска в атмосфере",
    )
    plt.plot(
        R * np.cos(angle) * 10 ** (-3),
        R * np.sin(angle) * 10 ** (-3),
        linewidth=2,
        color="black",
        label="Поверхность Титана",
    )
    plt.plot(
        (p0 / (1 + e0 * np.cos(psi_sh)) * np.cos(psi_sh)) * 10 ** (-3),
        p0 / (1 + e0 * np.cos(psi_sh)) * np.sin(psi_sh) * 10 ** (-3),
        marker="o",
        color="black",
        markersize=5,
        linestyle="",
        label="Точка схода СА с опорной орбиты",
    )
    plt.plot(
        r_vh * np.cos(Yfi1[0]) * 10 ** (-3),
        r_vh * np.sin(Yfi1[0]) * 10 ** (-3),
        marker=".",
        color="black",
        markersize=7,
        linestyle="",
        label="Точка входа СА в атмосферу",
    )
    plt.title(title)
    plt.legend(loc="best", prop={"size": 6})
    plt.xlabel("км"), plt.ylabel("км")
    plt.savefig("Схема полета")


# TRAECTORIYA_SPUSKA
def flight_trajectory_fig(title, titan_lw, marksize):
    """Podgotovka figuri dlya risovaniya kuska trajectorii spuska v atmosphere

    Parameters
    ----------
    title : string
        nazvanie grafika
    titan_lw : int
        tolshina linii grafika
    marksize : int
        razmer tochek na grafike
    """
    fig1, ax1 = plt.subplots()
    fig1.set_facecolor("white")
    plt.axis("equal")
    angle = np.linspace(0, 2 * np.pi, 1000)
    plt.plot(
        R * np.cos(angle) * 10 ** (-3),
        R * np.sin(angle) * 10 ** (-3),
        linewidth=titan_lw,
        color="black",
        linestyle="dashed",
        label="Поверхность Титана",
    )
    plt.plot(
        (R + h_vh) * np.cos(angle) * 10 ** (-3),
        (R + h_vh) * np.sin(angle) * 10 ** (-3),
        linestyle=":",
        linewidth=titan_lw / 2,
        color="black",
        label="Граница атмосферы",
    )
    plt.plot(
        r_vh * np.cos(psi_vh) * 10 ** (-3),
        r_vh * np.sin(psi_vh) * 10 ** (-3),
        marker="o",
        color="black",
        markersize=marksize,
        linestyle="",
        label="Точка входа СА в атмосферу",
    )
    plt.title(title)
    plt.xlabel("км"), plt.ylabel("км")


def flight_trajectory(Yh, Yfi, label, traj_lw, ls, marksize, marklabel):
    """Risovanie kuska trajectorii spuska v atmosphere

    Parameters
    ----------
    Yh : list of float
        visota poleta
    Yfi : list of float
        ugol poleta
    label : string
        parametr grafika dlya matplotlib
    traj_lw : int
        parametr grafika dlya matplotlib
    ls : string
        parametr grafika dlya matplotlib
    marksize : int
        parametr grafika dlya matplotlib
    marklabel : string
        parametr grafika dlya matplotlib
    """
    plt.plot(
        (Yh + R) * np.cos(Yfi) * 10 ** (-3),
        (Yh + R) * np.sin(Yfi) * 10 ** (-3),
        linestyle=ls,
        color=clr[num],
        linewidth=traj_lw,
        label=label,
    )
    plt.plot(
        R * np.cos(Yfi[-1]) * 10 ** (-3),
        R * np.sin(Yfi[-1]) * 10 ** (-3),
        marker="X",
        markersize=marksize,
        color="black",
        linestyle="",
        label=marklabel,
    )

    if (Yh[-1] + R) * np.cos(Yfi[-1]) > (h_vh + R) * np.cos(Yfi[0]):
        Xlim1 = (h_vh + R) * np.cos(Yfi[0]) * 10 ** (-3) - 150
        Xlim2 = (Yh[-1] + R) * np.cos(Yfi[-1]) * 10 ** (-3) + 150
    else:
        Xlim1 = (Yh[-1] + R) * np.cos(Yfi[-1]) * 10 ** (-3) - 150
        Xlim2 = (h_vh + R) * np.cos(Yfi[0]) * 10 ** (-3) + 150

    if (Yh[-1] + R) * np.sin(Yfi[-1]) > (h_vh + R) * np.sin(Yfi[0]):
        Ylim1 = (h_vh + R) * np.sin(Yfi[0]) * 10 ** (-3) - 150
        Ylim2 = (Yh[-1] + R) * np.sin(Yfi[-1]) * 10 ** (-3) + 150
    else:
        Ylim1 = (Yh[-1] + R) * np.sin(Yfi[-1]) * 10 ** (-3) - 150
        Ylim2 = (h_vh + R) * np.sin(Yfi[0]) * 10 ** (-3) + 150
    plt.xlim(Xlim1, Xlim2)
    plt.ylim(Ylim1, Ylim2)
    plt.savefig("Траектория спуска")


# SCOROST
def velocity_time_fig(title):
    """Podgotovka figuri dlya grafika scorosti

    Parameters
    ----------
    title : string
        nazvanie grafika
    """
    fig1, ax1 = plt.subplots()
    fig1.set_facecolor("white")
    plt.title(title)
    plt.xlabel("Время, с"), plt.ylabel("Скорость, м/с")


def velocity_time(T, Yv, ls, lw, label):
    """Risovanie grafika scorosti

    Parameters
    ----------
    T : np.array
        massiv vremeni
    Yv : np.array
        massiv skorostey
    ls : string
        parametr dlya matplotlib
    lw : int
        parametr dlya matplotlib
    label : string
        dop string dlya nazvaniya (s uchetom/ bez ucheta unosa TZP)
    """
    plt.plot(T, Yv, linewidth=lw, linestyle=ls, color=clr[num], label=label)
    plt.xlim(0, T[-1])
    plt.savefig("Скорость")


# UGOL_NAKLONA_VECT_SKOROSTI
def veloangle_time_fig(title, fs):
    """Podgotovka figuri dlya grafika ugla naklona scorosti

    Parameters
    ----------
    title : string
        nazvanie grafika
    fs : int
        razmer shrifta (currently unused)
    """
    fig1, ax1 = plt.subplots()
    fig1.set_facecolor("white")
    plt.title(title, fontsize=fs)
    plt.xlabel("Время, с"), plt.ylabel("Угол, град")


def veloangle_time(T, Ytheta, ls, lw, label):
    """Risovanie grafika ugla naklona skorosti

    Parameters
    ----------
    T : np.array
        massiv vremeni
    Ytheta : np.array
        massiv uglov naklona skorostey
    ls : string
        parametr dlya matplotlib
    lw : int
        parametr dlya matplotlib
    label : string
        dop string dlya nazvaniya (s uchetom/ bez ucheta unosa TZP)
    """
    plt.plot(
        T, np.degrees(Ytheta), linestyle=ls, linewidth=lw, color=clr[num], label=label
    )
    plt.xlim(0, T[-1])
    plt.savefig("Угол наклона вектора скорости")


# VISOTA_POLETA
def altitude_time_fig(title):
    """Podgotovka figuri dlya grafika visoti poleta

    Parameters
    ----------
    title : string
        nazvanie
    """
    fig1, ax1 = plt.subplots()
    fig1.set_facecolor("white")
    plt.title(title)
    plt.xlabel("Время, с"), plt.ylabel("Высота, км")


def altitude_time(T, Yh, ls, lw, label):
    """Risovanie grafika ugla naklona skorosti

    Parameters
    ----------
    T : np.array
        massiv vremeni
    Yh : np.array
        massiv visoti poleta
    ls : string
        parametr dlya matplotlib
    lw : int
        parametr dlya matplotlib
    label : string
        dop string dlya nazvaniya (s uchetom/ bez ucheta unosa TZP)
    """
    plt.plot(
        T, Yh * 10 ** (-3), linewidth=lw, linestyle=ls, color=clr[num], label=label
    )
    plt.xlim(0, T[-1])
    plt.savefig("Высота полета")


# DALNOST_POLETA
def length_time_fig(title):
    """Podgotovka figuri dlya grafika dalnosti poleta

    Parameters
    ----------
    title : string
        nazvanie
    """
    fig1, ax1 = plt.subplots()
    fig1.set_facecolor("white")
    plt.title(title)
    plt.xlabel("Время, с"), plt.ylabel("Дальность, км")


def length_time(T, YL, ls, lw, label):
    """Risovanie grafika ugla naklona skorosti

    Parameters
    ----------
    T : np.array
        massiv vremeni
    YL : np.array
        massiv dalnosti poleta
    ls : string
        parametr dlya matplotlib
    lw : int
        parametr dlya matplotlib
    label : string
        dop string dlya nazvaniya (s uchetom/ bez ucheta unosa TZP)
    """
    plt.plot(
        T, YL * 10 ** (-3), linewidth=lw, linestyle=ls, color=clr[num], label=label
    )
    plt.xlim(0, T[-1])
    plt.savefig("Дальность полета")


# PEREGRUZKA
def n_time_fig(title):
    """Podgotovka figuri dlya grafika peregruzki

    Parameters
    ----------
    title : string
        nazvanie
    """
    fig1, ax1 = plt.subplots()
    fig1.set_facecolor("white")
    plt.title(title)
    plt.xlabel("Время, с"), plt.ylabel("Перегрузка")


def n_time(T, Yn, ls, lw, dot, marksize, label):
    """Risovanie grafika ugla naklona skorosti

    Parameters
    ----------
    T : np.array
        massiv vremeni
    Yn : np.array
        massiv peregruzki
    ls : string
        parametr dlya matplotlib
    lw : int
        parametr dlya matplotlib
    dot : int
        perekluchatel
    marksize : int
        parametr dlya matplotlib
    label : string
        dop string dlya nazvaniya (s uchetom/ bez ucheta unosa TZP)
    """
    plt.plot(T, Yn, linewidth=lw, linestyle=ls, color=clr[num], label=label)
    if dot != 0:
        plt.plot(
            T[imax],
            Yn[imax],
            marker="o",
            color="black",
            markersize=marksize,
            linestyle="",
            label="nmax = %.2f , " "t = %.0f с" % (Yn[imax], T[imax]),
        )
    else:
        1
    plt.xlim(0, T[-1])
    plt.savefig("Перегрузка")


# DEQSTV_SILI
def force_sum_fig(title):
    """Podgotovka figuri dlya grafika sil

    Parameters
    ----------
    title : string
        nazvanie
    """
    fig1, ax1 = plt.subplots()
    fig1.set_facecolor("white")
    plt.title(title)
    plt.xlabel("Время, с"), plt.ylabel("Сила, кН")


def force_sum(T, Yh, Yv, Ytheta, lw, style):
    """Risovanie grafika sil

    Parameters
    ----------
    T : np.array
        massiv vremeni
    Yh : np.array
        massiv visoti
    Yv : np.array
        massiv skorosti
    Ytheta : np.array
        massiv ugla naklona skorosti
    lw : int
        parametr dlya matplotlib
    style : int
        parametr dlya matplotlib
    """
    if style == 1:
        label = ["Без учета уноса ТЗП", "", "", "", ""]
        ls = "dashed"
        cl = ["dimgrey", "dimgrey", "dimgrey", "dimgrey", "dimgrey"]
    else:
        label = [
            "Сила лобового сопротивления X",
            "Сила тяжести Gx",
            "Подъемная сила Y",
            "Сила тяжести Gy",
            "Центробежная сила",
        ]
        ls = "solid"
        cl = ["tab:blue", "tab:orange", "tab:green", "tab:red", "tab:purple"]
    plt.plot(
        T,
        -Cx * Sm * spline(Yh) * Yv**2 * 10 ** (-3) / 2,
        linewidth=lw,
        color=cl[0],
        linestyle=ls,
        label=label[0],
    )
    plt.plot(
        T,
        -M * f * np.sin(Ytheta) * 10 ** (-3) / (Yh + R) ** 2,
        linewidth=lw,
        color=cl[1],
        linestyle=ls,
        label=label[1],
    )
    plt.plot(
        T,
        Cy * Sm * spline(Yh) * Yv**2 * 10 ** (-3) / 2,
        linewidth=lw,
        color=cl[2],
        linestyle=ls,
        label=label[2],
    )
    plt.plot(
        T,
        -M * f * np.cos(Ytheta) * 10 ** (-3) / (Yh + R) ** 2,
        linewidth=lw,
        color=cl[3],
        linestyle=ls,
        label=label[3],
    )
    plt.plot(
        T,
        M * Yv**2 * np.cos(Ytheta) * 10 ** (-3) / (Yh + R),
        linewidth=lw,
        color=cl[4],
        linestyle=ls,
        label=label[4],
    )
    plt.xlim(0, T[-1])
    plt.savefig("Силы действующие на СА")


# UNOSIMAYA_TZP
def mass_time_fig(title):
    """Podgotovka figuri dlya grafika unosa TZP

    Parameters
    ----------
    title : string
        nazvanie
    """
    fig1, ax1 = plt.subplots()
    fig1.set_facecolor("white")
    plt.title(title)
    plt.xlabel("Время, с"), plt.ylabel("Масса, кг")


def mass_time(T, Ym):
    """Risovanie grafika unosa TZP

    Parameters
    ----------
    T : np.array
        massiv vremeni
    Ym : np.array
        massiv massi sgorevshey TZP
    """
    plt.plot(T, Ym * Sm, color="black", linewidth=1.5)
    plt.xlim(0, T[-1])
    plt.savefig("Уносимая ТЗП")


# TEMPERATURA
def temp_time_fig(title):
    """Podgotovka figuri dlya grafika temperaturi

    Parameters
    ----------
    title : string
        nazvanie
    """
    fig1, ax1 = plt.subplots()
    fig1.set_facecolor("white")
    plt.title(title)
    plt.xlabel("Время, с"), plt.ylabel("Температура, К")


def temp_time(T, Yh, Yv, k1, k2, kmax, lw, marksize, label):
    """Risovanie grafika temperaturi

    Parameters
    ----------
    T : np.array
        massiv vremeni
    Yh : np.array
        massiv visoti
    Yv : np.array
        massiv skorosti
    k1 : int
        Temp nachala unosa TZP
    k2 : int
        Temp konec unosa TZP
    kmax : int
        Max temp
    lw : int
        parametr dlya matplotlib
    marksize : int
        parametr dlya matplotlib
    label : string
        dop string dlya nazvaniya
    """
    if k1 == 0 or k2 == 0:
        a = False
    else:
        a = True
    if label == "":
        label = "Температура ТЗП"
    plt.plot(T, Tp(Yh, Yv), linewidth=lw, label=label, color=clr[num])
    plt.plot(
        T[kmax],
        Tp(Yh[kmax], Yv[kmax]),
        marker="o",
        color="black",
        markersize=marksize,
        linestyle="",
        label="Tmax = %.0f К, " "t = %.0f с" % (Tp(Yh[kmax], Yv[kmax]), T[kmax]),
    )
    if a:
        plt.plot(
            T[k1],
            Tp(Yh[k1], Yv[k1]),
            marker="o",
            color="black",
            markersize=marksize,
            linestyle="",
            label="Начало разрушения ТЗП - %.0f с" % T[k1],
        )
        plt.plot(
            T[k2],
            Tp(Yh[k2], Yv[k2]),
            marker="o",
            color="black",
            markersize=marksize,
            linestyle="",
            label="Конец разрушения ТЗП - %.0f с" % T[k2],
        )
        plt.plot(
            T,
            np.linspace(Td_tz, Td_tz, len(T)),
            linewidth=0.75,
            linestyle="dashed",
            color="black",
            label="Температура " "разрушения ТЗП",
        )
        plt.axvline(T[k1], color="black", linewidth=0.6, linestyle="dashdot")
        plt.axvline(T[k2], color="black", linewidth=0.6, linestyle="dashdot")
    plt.xlim(0, T[-1])
    plt.savefig("Температура")


# KONVEKCIYA
def heatk_time_fig(title):
    """Podgotovka figuri dlya grafika konvectivnogo teploobmena

    Parameters
    ----------
    title : string
        nazvanie
    """
    fig1, ax1 = plt.subplots()
    fig1.set_facecolor("white")
    plt.title(title)
    plt.xlabel("Время, с"), plt.ylabel("Плотность потока, МДж/(с*м^2)")


def heatk_time(T, Yh, Yv):
    """Risovanie grafika konvectivnogo teploobmena

    Parameters
    ----------
    T : np.array
        massiv vremeni
    Yh : np.array
        massiv visoti
    Yv : np.array
        massiv skorosti
    """
    plt.plot(T, qk(Yh, Yv) * 10 ** (-6), color="black", linewidth=1.5)
    plt.xlim(0, T[-1])
    plt.savefig("Конвекция")


# RADIACIYA
def heatr_time_fig(title):
    """Podgotovka figuri dlya grafika radiacionnogo teploobmena

    Parameters
    ----------
    title : string
        nazvanie
    """
    fig1, ax1 = plt.subplots()
    fig1.set_facecolor("white")
    plt.title(title)
    plt.xlabel("Время, с"), plt.ylabel("Плотность потока, мкДж/(с*м^2)")


def heatr_time(T, Yh, Yv):
    """Risovanie grafika radiacionnogo teploobmena

    Parameters
    ----------
    T : np.array
        massiv vremeni
    Yh : np.array
        massiv visoti
    Yv : np.array
        massiv skorosti
    """
    plt.plot(T, qr(Yh, Yv) * 10**6, color="black", linewidth=1.5)
    plt.xlim(0, T[-1])
    plt.savefig("Радиация")


# Resheniye systemy uravneniya

ballistic(v_vh, theta_vh, psi_vh)

y0, t0 = [v_vh, theta_vh, 0, h_vh], 0
tmax = 100000
n_steps = 10**5
ODE = ode(equations1)
ODE.set_integrator("dopri5", nsteps=n_steps)
ODE.set_solout(step1)
ODE.set_initial_value(y0, t0)
ODE.integrate(tmax)

Yv1, Ytheta1, YL1, Yh1 = Y1[:, 0], Y1[:, 1], Y1[:, 2], Y1[:, 3]
Yfi1 = Yfi_find(T1, psi_vh, Yv1, Yh1)
Yn1 = Yn_nmax(Cx, Sm, Yh1, Yv1, K)

clr = ["black", "red", "orange", "green", "blue", "limegreen", "gold", "magenta"]
num = 0

if os.path.isdir("dz/titan_lander_50/graph"):
    shutil.rmtree("dz/titan_lander_50/graph")
os.mkdir("dz/titan_lander_50/graph")
os.chdir("dz/titan_lander_50/graph")


flight_scheme("Траектория маневра без учета уноса ТЗП", Yh1, Yfi1)
flight_trajectory_fig("Траектория спуска СА без учета уноса ТЗП", 1, 5)
flight_trajectory(Yh1, Yfi1, "", 1, "solid", 5, "Точка приземления СА")
plt.legend(loc="best", prop={"size": 8})

velocity_time_fig(
    "Зависимость величины вектора скорости от времени\n" "без учета уноса ТЗП"
)
plt.grid(color="gray", linestyle="dotted")
velocity_time(T1, Yv1, "solid", 1.5, ""), plt.ylim(ymin=0)

veloangle_time_fig(
    "Зависимость угла наклона вектора скорости к местному\n"
    "горизонту от времени без учета уноса ТЗП",
    None,
)
plt.grid(color="gray", linestyle="dotted")
veloangle_time(T1, Ytheta1, "solid", 1.5, "")

altitude_time_fig("Зависимость высоты полета от времени\n" "без учета уноса ТЗП")
plt.grid(color="gray", linestyle="dotted")
altitude_time(T1, Yh1, "solid", 1.5, ""), plt.ylim(ymin=0)

length_time_fig("Зависимость дальности полета от времени\n" "без учета уноса ТЗП")
plt.grid(color="gray", linestyle="dotted")
length_time(T1, YL1, "solid", 1.5, ""), plt.ylim(ymin=0)

n_time_fig("Зависимость коэффициента перегрузки от времени\n" "без учета уноса ТЗП")
plt.grid(color="gray", linestyle="dotted")
n_time(T1, Yn1, "solid", 1.5, 0, 5, ""), plt.ylim(ymin=0)

os.chdir("../../../")


"РАСЧЕТ И ВЫВОД СЛУЧАЯ С УЧЕТОМ УНОСА ТЗП"


if os.path.isdir("dz/titan_lander_50/graph_u"):
    shutil.rmtree("dz/titan_lander_50/graph_u")
os.mkdir("dz/titan_lander_50/graph_u")
os.chdir("dz/titan_lander_50/graph_u")

y0, t0 = [v_vh, theta_vh, 0, h_vh, 0], 0
ODE = ode(equations2)
ODE.set_integrator("dopri5", nsteps=n_steps)
ODE.set_solout(step2)
ODE.set_initial_value(y0, t0)
ODE.integrate(tmax)

Yv2, Ytheta2, YL2, Yh2, Ym2 = Y2[:, 0], Y2[:, 1], Y2[:, 2], Y2[:, 3], Y2[:, 4]
temp_points(T2, Yh2, Yv2)
Yfi2 = Yfi_find(T2, psi_vh, Yv2, Yh2)
Yn2 = Yn_nmax(Cx, Sm, Yh2, Yv2, K)


flight_scheme("Траектория маневра с учетом уноса ТЗП", Yh2, Yfi2)

flight_trajectory_fig("Траектория спуска СА с учетом уноса ТЗП", 1, 5)
flight_trajectory(Yh2, Yfi2, "", 1, "solid", 5, "Точка приземления СА")
plt.legend(loc="best", prop={"size": 8})

velocity_time_fig(
    "Зависимость величины вектора скорости от времени\n" "с учетом уноса ТЗП"
)
plt.grid(color="gray", linestyle="dotted")
velocity_time(T2, Yv2, "solid", 1.5, ""), plt.ylim(ymin=0)

veloangle_time_fig(
    "Зависимость угла наклона вектора скорости к местному\n"
    "горизонту от времени с учетом уноса ТЗП",
    None,
)
plt.grid(color="gray", linestyle="dotted")
veloangle_time(T2, Ytheta2, "solid", 1.5, "")

altitude_time_fig("Зависимость высоты полета от времени\n" "с учетом уноса ТЗП")
plt.grid(color="gray", linestyle="dotted")
altitude_time(T2, Yh2, "solid", 1.5, ""), plt.ylim(ymin=0)

length_time_fig("Зависимость дальности полета от времени\n" "с учетом уноса ТЗП")
plt.grid(color="gray", linestyle="dotted")
length_time(T2, YL2, "solid", 1.5, ""), plt.ylim(ymin=0)

n_time_fig("Зависимость коэффициента перегрузки от времени\n" "с учетом уноса ТЗП")
plt.grid(color="gray", linestyle="dotted")
n_time(T2, Yn2, "solid", 1.5, 0, 5, ""), plt.ylim(ymin=0)

mass_time_fig("Зависимость уносимой массы ТЗП от времени")
plt.grid(color="gray", linestyle="dotted")
mass_time(T2, Ym2), plt.ylim(ymin=0)

temp_time_fig("Зависимость температуры ТЗП от времени")
plt.grid(color="gray", linestyle="dotted")
temp_time(T2, Yh2, Yv2, k1, k2, kmax, 1, 5, "")
plt.legend(loc="best", prop={"size": 8}), plt.ylim(ymin=0)

heatk_time_fig("Зависимость плотности конвективного\nтеплового потока от времени")
plt.grid(color="gray", linestyle="dotted")
heatk_time(T2, Yh2, Yv2), plt.ylim(ymin=0)

heatr_time_fig("Зависимость плотности радиационного\nтеплового потока от времени")
plt.grid(color="gray", linestyle="dotted")
heatr_time(T2, Yh2, Yv2), plt.ylim(ymin=0)

os.chdir("../../../")


if os.path.isdir("dz/titan_lander_50/graph_sravn"):
    shutil.rmtree("dz/titan_lander_50/graph_sravn")
os.mkdir("dz/titan_lander_50/graph_sravn")
os.chdir("dz/titan_lander_50/graph_sravn")


flight_trajectory_fig("Траектории спуска СА", 2, 6)
flight_trajectory(Yh1, Yfi1, "Без учета уноса ТЗП", 0.75, "solid", 6, "")
flight_trajectory(Yh2, Yfi2, "С учетом уноса ТЗП", 0.75, "dashed", 6, "")
plt.legend(loc="best", prop={"size": 8})

velocity_time_fig("Зависимости величины вектора скорости от времени")
velocity_time(T1, Yv1, "solid", 0.5, "Без учета уноса ТЗП")
velocity_time(T2, Yv2, "dashed", 0.5, "С учетом уноса ТЗП")
plt.legend(loc="best", prop={"size": 8}), plt.ylim(ymin=0)

veloangle_time_fig(
    "Зависимости угла наклона вектора скорости к местному\n" "горизонту от времени",
    None,
)
veloangle_time(T1, Ytheta1, "solid", 0.5, "Без учета уноса ТЗП")
veloangle_time(T2, Ytheta2, "dashed", 0.5, "С учетом уноса ТЗП")
plt.legend(loc="best", prop={"size": 8})

altitude_time_fig("Зависимости высоты полета от времени")
altitude_time(T1, Yh1, "solid", 0.5, "Без учета уноса ТЗП")
altitude_time(T2, Yh2, "dashed", 0.5, "С учетом уноса ТЗП")
plt.legend(loc="best", prop={"size": 8}), plt.ylim(ymin=0)

length_time_fig("Зависимости дальности полета от времени")
length_time(T1, YL1, "solid", 0.5, "Без учета уноса ТЗП")
length_time(T2, YL2, "dashed", 0.5, "С учетом уноса ТЗП")
plt.legend(loc="best", prop={"size": 8}), plt.ylim(ymin=0)

n_time_fig("Зависимости коэффициента перегрузки от времени")
Yn_nmax(Cx, Sm, Yh1, Yv1, K)
n_time(T1, Yn1, "solid", 0.5, 1, 3, "Без учета уноса ТЗП")
Yn_nmax(Cx, Sm, Yh2, Yv2, K)
n_time(T2, Yn2, "dashed", 0.5, 1, 3, "С учетом уноса ТЗП")
plt.legend(loc="best", prop={"size": 8}), plt.ylim(ymin=0)


force_sum_fig("Динамика действующих на СА сил")
# force_sum(T1, Yh1, Yv1, Ytheta1, 0.5, 1)
force_sum(T2, Yh2, Yv2, Ytheta2, 0.5, 2)
plt.legend(loc="best", prop={"size": 6})

os.chdir("../")

doc = DocxTemplate("Report.docx")
pict1 = InlineImage(doc, "graph/Высота полета.png", Cm(12))
pict2 = InlineImage(doc, "graph/Дальность полета.png", Cm(12))
pict3 = InlineImage(doc, "graph/Перегрузка.png", Cm(12))
pict4 = InlineImage(doc, "graph/Скорость.png", Cm(12))
pict5 = InlineImage(doc, "graph/Схема полета.png", Cm(12))
pict6 = InlineImage(doc, "graph/Траектория спуска.png", Cm(12))
pict7 = InlineImage(doc, "graph/Угол наклона вектора скорости.png", Cm(12))
pict8 = InlineImage(doc, "graph_u/Высота полета.png", Cm(12))
pict9 = InlineImage(doc, "graph_u/Дальность полета.png", Cm(12))
pict10 = InlineImage(doc, "graph_u/Перегрузка.png", Cm(12))
pict11 = InlineImage(doc, "graph_u/Скорость.png", Cm(12))
pict12 = InlineImage(doc, "graph_u/Схема полета.png", Cm(12))
pict13 = InlineImage(doc, "graph_u/Траектория спуска.png", Cm(12))
pict14 = InlineImage(doc, "graph_u/Угол наклона вектора скорости.png", Cm(12))
pict15 = InlineImage(doc, "graph_u/Конвекция.png", Cm(12))
pict16 = InlineImage(doc, "graph_u/Радиация.png", Cm(12))
pict17 = InlineImage(doc, "graph_u/Температура.png", Cm(12))
pict18 = InlineImage(doc, "graph_u/Уносимая ТЗП.png", Cm(12))
context = {
    "Vis1": Yv1[-1],
    "Vis2": Yv2[-1],
    "Dal1": YL1[-1] * 10 ** (-3),
    "Dal2": YL2[-1] * 10 ** (-3),
    "n1": Yn1[imax],
    "n2": Yn2[imax],
    "Pict1": pict1,
    "Pict2": pict2,
    "Pict3": pict3,
    "Pict4": pict4,
    "Pict5": pict5,
    "Pict6": pict6,
    "Pict7": pict7,
    "Pict8": pict8,
    "Pict9": pict9,
    "Pict10": pict10,
    "Pict11": pict11,
    "Pict12": pict12,
    "Pict13": pict13,
    "Pict14": pict14,
    "Pict15": pict15,
    "Pict16": pict16,
    "Pict17": pict17,
    "Pict18": pict18,
}
doc.render(context)
doc.save("Report_done.docx")
